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ABSTRACT 


Investigation  of  steady  and  unsteady  flowfields  over  airfoils  is  an  active 
area  of  current  computational  and  experimental  research.  In  this  study,  the 
compressible,  viscous  flow  over  a  single  and  multi-element  airfoil  is  numerically 
simulated  by  solving  the  Navier  Stokes  equations.  The  motivation  for  this  work 
includes  interest  in  studying  the  effects  of  a  stationaiy/flapping  airfoil  combination 
in  tandem  configuration.  A  smgle-block  Navier-  Stokes  (NS)  solver  is  employed 
to  compute  unsteady  flowfields.  Turbulence  is  treated  using  the  Baldwin-Lomax 
turbulence  model.  A  single  C-grid  is  generated  and  it  is  partially  distorted  to 
simulate  the  flapping  motion.  Numerical  solutions  are  obtained  for  flows  at  a  fixed 
angle  of  attack  and  for  unsteady  flows  over  flapping  airfoils.  The  numerical 
solutions  agree  well  with  the  experimental  data.  The  difficulties  faced  dining  the 
study  are  discussed  and  future  improvements  are  suggested. 
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INTRODUCTION 


Studies  of  steady  and  unsteady  flowfields  over  multi-element  airfoils  is  an  active 
area  of  numerical  and  experimental  research.  While  notable  progress  is  being  made  in  the 
development  of  computing  methods  to  analyze  flows  around  single  airfoils,  methods  to 
compute  more  complex  geometric  configurations  are  now  being  developed.  In  this  study, 
a  Navier  Stokes  solver  was  modified  to  accommodate  the  complexity  of  a  multi-element 
configuration  to  study  the  thrust  generating  effect  of  a  stationary  /flapping  airfoil 
combination. 

Thrust  generation  due  to  airfoil  flapping  was  recognized  by  early  researchers,  such 
as  Knoller  [Ref  1],  in  explaining  the  bird’s  ability  to  generate  a  propulsive  force  by  means 
of  flapping  their  wing.  Previous  experimental  investigations  of  unsteady  flows  over 
oscillating  airfoils  by  Neace  [Ref  2]  ,Tuncer  and  Platzer  [Ref  3],  and  Dohring  [Ref  4] 
have  shown  a  propulsive  effect  and  the  associated  eflBciencies.  Figure  1  shows  how  the 
airfoil  motion  creates  an  induced  velocity  component.  This  allows  for  the  generation  of  a 
forward  or  thrust  component  of  force. 


In  the  process  of  this  research,  a  single  block  Navier  Stokes  solver  was  first 
validated  on  a  Ames-01  airfoil.  The  solver  was  subsequently  modified  to  handle  multi¬ 
element  airfoil  configurations  utilizing  a  single  structured  C-grid.  Computational  grids 
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were  generated  utilizing  GRIDGEN  .  Other  grid  generating  tools  were  evaluated,  but 
only  GRIDGEN  provided  the  flexibility  to  properly  distribute  the  gridlines  in  critical  areas 
The  initial  grid  configuration  utilizing  NACA  4412  and  4415  airfoils  within  a  wind  tunnel 
boundary  was  used  in  the  validation  process  of  the  code.  Two  other  configurations  were 
studied  under  unsteady  conditions  .  In  the  unsteady  flow  case  the  leading  airfoil  is  kept 
stationary  while  the  trailing  airfoil  undergoes  a  flapping  motion.  In  Case  I,  flap  is  placed 
at  .5c  aft  of  the  main  airfoil  and  in  Case  II  at  .25c  aft  of  the  main  airfoil  These  two  cases 
were  studied  at  different  flapping  frequencies  and  amplitudes. 
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n.  METHOD 


A.  GRH)  GENERATION 

1.  Code  Validation  Grid 

During  the  grid  generation  process  proper  distribution  of  grid  points  around  the 
airfoil’s  leading  edge  and  trailing  edge,  orthogonality  of  the  grid  lines  on  the  surface,  and 
smooth  variations  of  the  grid  density  are  the  most  critical  factors.  In  order  to  accurately 
calculate  the  flow  gradients  in  the  direction  normal  to  the  surface  in  the  boundary  layer,  it 
is  necessary  to  make  the  normal  grid  spacing  very  fine  close  to  the  sohd  surfaces.  Single 
C-type  grids  for  NACA  4412  (Main  airfoil)  and  NACA  4415  (Flap)  were  generated  for 
the  code  validation  process  modeling  the  configuration  in  Figure  2.  The  grids  included  in 
this  report  were  generated  using  GRIDGEN  software  packages. 

Adair  and  Home  [Ref  5]  conducted  experiments  with  this  configuration  in  the  7 
by  10  Foot  Wind  Tunnel  at  NASA  Ames  Research  Center.  In  the  experiment,  the  chord 
length  of  the  main  airfoil  is  .90  m  and  that  of  the  flap  is  .36  m.  The  geometric  location  of 
the  flap  relative  to  the  main  airfoil  was  described  by  the  flap  gap  (FG) ,  the  flap  overlap 
(FO)  and  the  flap  deflection  (5f).  For  this  case  FG  =  0.035c ,  FO  =  0.028c ,  6f  =  21.8° 
and  the  main  airfoil  angle  of  attack  (a)  was  set  to  8.2°. 


Figure  2  Experimental  Configuration  Adair  and  Home  [Ref.  5] 
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The  proximity  of  the  airfoils,  and  the  sensitivity  of  the  solver  puts  a  heavy 
constraint  on  the  computational  grid  around  the  downstream  airfoil.  In  this  case  it  was 
not  possible  to  duplicate  the  exact  positioning  of  the  flap,  Instead,  the  flap  was  positioned 
at  0.015c  aft  of  the  main  airfoil,  as  shown  in  Figure  3,  and  a  parametric  study  on  the  flap 
overlap  and  gap  was  conducted  to  verify  proper  trends  and  extrapolate  verification.  To 
improve  the  resolution  and  distribution  of  the  gridlines  in  the  wake  of 


the  main  airfoil  GRIDGEN’s  elliptical  solver  was  used.  One  of  the  main  reasons  for  using 
an  elliptic  solver  is  to  force  the  grid  to  have  a  smooth  variation  while  maintaining  the  grid 
as  body  fitted.  The  other  reason  is  the  ability  of  the  solution  to  provide  clustering  and 
orthogonality  to  the  grid  based  on  the  control  function  selected.  GRIDGEN  offers  a 
choice  of  four  control  function  types  that  influences  the  grid  distribution  and 
orthogonality.  Among  those,  the  Laplace  control  function  with  a  ,7  relaxation  parameter 
was  the  standard  selection  for  this  research.  This  type  of  function  is  the  least 
computationally  demanding  and  the  most  stable.  This  method  produces  the  smoothest 
possible  grid.  For  a  detailed  description  of  all  related  parameters  see  the  GRIDGEN 
Manual  [Ref  6].  See  sample  of  GRIDGEN  input  file  and  elliptical  solver  parameters  in 
the  Appendix  B. 
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The  following  is  a  list  of  parameters  and  procedures  used  to  run  the  elliptical  solver: 

•  Relaxation  factor  =  0.7 

•  Laplace  control  function 

•  Set  EC’s  on  the  airfoil  surface  to  orthogonality.  By  default  all  points  on  the 
edges  will  remain  fixed  as  the  elliptical  solver  is  run. 

•  Turn  foreground  EC’s  on  (elliptical  solver  menu)  with  default  values. 

•  Run  the  solver  for  critical  areas  separately,  this  will  expedite  the  process  and 
will  result  in  a  smoother  grid. 

•  Run  the  solver  until  an  acceptable  grid  is  achieved. 

Figure  4  shows  the  grid  distribution  in  the  wake  area  for  the  initial  configuration 
and  the  up-flap  configuration.  For  the  aft-flap  configuration  the  flap  was  moved  aft  0.03c, 
for  the  up-flap  configuration  the  flap  was  moved  up  from  the  aft  position  0.03c.  Due  to 
the  positioning  of  the  flap  in  the  wake  of  the  main  airfoil,  (up-flap  configuration)  the 
density  of  the  grid  lines  between  airfoils  was  increased.  The  grid  dimensions  for  the  initial 
and  up-flap  configurations  were  361  x  100  and  463  x  100,  respectively. 


Figure  5  shows  the  full  wind  tunnel  model  for  the  aft-flap  configuration  used 
during  the  validation  process. 
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Figure  5  Aft-flap  Wind  Tunnel  Grid 


2.  Unsteady  Case  I  and  11  Grids 

For  unsteady  computations  two  separate  grids  were  generated  as  Case  I  and  II. 
NACA  4412  and  4415  airfoils  were  used  for  these  computations.  The  airfoil  combination 
was  set  in  a  tandem  configuration  with  a  gap  between  airfoils  of  .5c  for  Case  I  and  .25c 
for  Case  II.  (Figure  6).  The  length  of  the  flap  (NACA  4415)  was  set  to  0.33c.  The  main 
airfoil  was  set  at  a  =  10"  and  the  flap  at  zero  angle  of  attack  with  respect  to  the  main 
airfoil  chord  line.  The  grid  dimensions  are  391  x  100  for  Case  I  and  371  x  100  for  Case 
II 

These  configurations  eliminate  some  of  the  complexities  of  the  initial  configuration 
grid  allowing  the  evaluation  of  the  oscillation  effects  without  the  influence  of  any  other 
variables.  In  addition,  having  the  trailing  airfoil  far  enough  from  the  main  airfoil  facilitates 
the  implementation  of  distortion  required  to  model  the  oscillation  process.  Details  of  the 
grid  distortion  are  discussed  in  the  next  section. 


Figure  6  Case  I  (top)  and  II  grids(bottom) 


B.  ALGORITHM 

1.  Navier-Stokes  Solver 

An  implicit,  thin  layer,  Navier-Stokes  solver  with  the  third  order  accurate  Osher's 
upwind  biased  flux  difference  splitting  scheme  is  employed  to  compute  the  flowfield  in  the 
computational  grid.  The  strong  conservation  law  form  of  the  2-D,  thin  layer  Navier- 
Stokes  equations  in  a  curvilinear  coordinate  system,  (^,5),  along  the  axial  and  normal 
direction  respectively,  is  given  as  follows 

46  +  4^  +  45  (2. 1) 

Where  Q  is  the  vector  of  conservative  variables,  l/J(  p,  pu,  pw,  e),  F  and  G  are  the 

inviscid  flux  vectors,  and  S  is  the  thin  layer  approximation  of  the  viscous  fluxes  in  the  q 
direction  normal  to  the  airfoil  surface: 
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where 


^  pU  ^ 

^  pW  ^ 

1 

puU  +  4;? 

r-i- 

puW+^^p 

J 

pwU  +  ^^p 

K{e  +  p)U-^pj 

l(e+pW-CpJ 

J  +  {ju  / 

+{iu  1 3)m^  + 


(2.2) 


^2  = 


Wj  =(u^  +  w^)/2+k:?t~^ 


(2.3) 


and  U.W axe  the  contravariant  velocity  components.  In  the  above  equations,  all 
dimensions  are  normalized  with  the  airfoil  chord  length,  c.  p  is  the  density  normalized 
with  the  free-stream  density  /?«;  u  and  w  are  the  Cartesian  velocity  components  in  the 
physical  domain,  which  are  normalized  with  the  free-stream  speed  of  sound  a^;  e  is  the 

total  energy  per  unit  volume  normalized  with  p^al^ ;  and  Pr  is  the  Prandtl  number.  The 
pressure  is  related  to  density  and  total  energy  through  the  equation  of  state  for  an  ideal 
%as,  p  =  {r-\)[e-p{u^ +w^)/2'\.  (2.4) 

The  flowfield  is  assumed  to  be  fully  turbulent.  The  turbulence  modeling  is  used  to  relate 
the  Reynolds  shear  stress  to  the  local  mean-velocity  gradient  allowing  numerical  flow  field 
calculation.  The  Baldwin—Lomax  algebraic  turbulence  model  is  currently  implemented. 
This  model  is  a  two  layer  eddy  viscosity  model  which  simulates  the  effect  of  turbulence  in 
terms  of  eddy  viscosity  coefficient.  A  complete  description  of  the  model  is  given  in 
Baldwin  [Ref  7]. 
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2.  Boundary  Conditions 

The  computational  domain  includes  both  airfoil  surfaces  and  extends  ten  chord 
lengths  away  from  the  airfoils.  Boundary  conditions  are  applied  on  the  airfoil  surfaces,  and 
at  the  farfield  boundaries.  On  the  airfoil  surfaces  the  no  slip  boundary  condition  is 
applied.  For  the  flapping  airfoil,  the  surface  fluid  velocity  is  set  equal  to  the  prescribed 
airfoil  flapping  velocity  so  that  the  no  shp  condition  is  satisfied.  Since  the  formulation  of 
the  Navier-Stokes  solver  is  based  on  an  inertial  frame  of  reference,  the  flapping  motion  of 
the  airfoil  is  implemented  by  moving  the  airfoil  and  the  computational  grid  around  it  in  the 
transverse  direction  as  described  by  the  frequency  and  the  amplitude  of  the  flapping 
motion  by: 

h  =  -Acos{cot) 

Where  A  denotes  the  mean  amplitude  normalized  with  the  chord  length,  o  is  the 
frequency  of  the  oscillatory  flapping  motion,  which  is  used  in  terms  of  reduced  frequency, 

k  =  (0cl2U^  (2-6) 

At  the  farfield  inflow  and  outflow  boundaries  the  flow  variables  are  evaluated 
using  the  zero  order  Riemann  invariant  extrapolation. 

For  the  code  validation  studies,  no  slip  boundary  condition  is  apphed  at  the  tunnel 

walls. 


3.  Unsteady-Motion  Routine 

The  features  pertinent  to  the  grid  distortion  in  the  NS  code  were  modified  so  that 
the  trailing  airfoil  can  move  in  the  cross-flow  direction  with  respect  to  the  leading  airfoil. 
During  the  flapping  motion  the  main  airfoil  remains  stationary  while  the  flap  oscillates  in 
the  cross  flow  direction.  This  motion  is  introduced  by  a  partial  grid  distortion.  The  wake 
region  between  airfoils  is  the  area  of  concern  since  the  distortion  is  gradually  introduced 
from  the  trailing  edge  of  the  main  airfoil  to  the  flap  leading  edge.  During  this  distortion 
the  integrity  of  the  grid  in  terms  of  orthogonality,  smoothness  and  clustering  of  grid  point 
should  remain  physically  valid  for  a  solution  to  converge.  In  order  to  maint^^in  a  smooth 
variation  of  grid  distortion  a  hyperbolic  tangent  is  used  as  the  transition  function  between 
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the  stationary  and  moving  part  of  the  grid..  The  grid’s  distortion  adjusts  smoothly  using 
this  function.  Figure  7  shows  this  function.  The  grid  motion  routine  (GRMOVE 
subroutine  of  the  NS  code)  is  given  in  the  Appendix  A. 


h  =  0.0 


Figure  7  Transition  Function-  Hyperbolic  Tangent 


m.  RESULTS  AND  DISCUSSION 


The  validity  of  the  code  was  first  tested  with  a  specified  configuration.  The  code 
was  tested  under  the  steady-state  conditions  described  below.  The  solution  was  then 
evaluated  for  a  number  of  configurations  and  compared  to  the  experimental  data. 
Following  the  validation  process  the  flowfield  with  the  flapping  trailing  airfoil  was  studied 
for  two  different  configurations. 

A.  CODE  VALIDATION 

1.  Conditions—Steady  State 

For  the  validation  process  the  flowfield  was  computed  at: 

•  a  =  8.2° 

•  Mach  =  .2 

•  Reynolds  number  1.8  x  10^,  fully  turbulent. 

The  solver  was  first  run  with  the  initial  wing-flap  configuration.  Convergence  was 
reached  at  1500  iterations  based  on  the  variations  in  the  flow  residuals  and  the 
aerodynamic  loads.  Similarly,  the  solver  was  run  for  the  aft-flap  and  up-aft  flap 
configurations.  A  typical  NS  input  file  is  shown  in  the  Appendix  A.  Figure  8  shows  the 
convergence  history  of  the  initial  case  and  Figure  9  shows  the  aerod5miamic  loads  for  the 
same  configuration. 


Figure  8  Convergence  History  Initial  Configuration 
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Iterations 


Figure  9  Aerodynamic  Loads  -  Initial  Configuration 

Figure  10  shows  the  pressure  coefficient  distribution  for  all  cases.  The  pressure  changes 
indicated  in  this  figure  for  the  different  cases  shows  a  general  trend  that  represents  a 


Figure  10  Computed  vs.  Experimental  Cp  Distribution 
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good  agreement  with  the  experimental  data.  Figure  1 1  summarizes  the  effects  of  the  gap 
on  the  separated  flap  flow.  This  figure  shows  how  the  flowfield  is  affected  as  the  flap  is 
lowered  or  moved  aft.  As  the  gap  is  increased  more  fluid  from  below  the  main  airfoil 
moves  into  the  flap  through  the  gap  pushing  the  flow  farther  away  from  the  flap  leading 
edge,  Alemdaroglu  [Ref  8].  As  the  gap  is  increased  the  flow  will  move  the  streamlines 
upward  affecting  the  flow  separation  over  the  flap.  This  condition  has  significant  effect  on 
the  pressure  distribution  over  the  aft-upper  surface  of  the  main  airfoil  as  noted  in  Figure 
10. 


Mass  Flux 


Figure  11  Computed  flowfield  -  Initial  Configuration  (top).  Aft-flap  at  8f=lS.  6° 

(bottom) 
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B.  UNSTEADY  CASE 


This  part  of  the  research  focused  on  the  aerodynamic  interaction  between  two 
NACA  airfoils  as  shown  in  Fi^re  6.  The  main  airfoil  (NACA  4412)  is  stationary  while 
the  flap  (NACA  441 5)  undergoes  a  flapping  motion.  The  flowfields  were  evaluated  under 
the  following  conditions: 

•  a  =  10  degrees 

•  Mach  =  .3 

•  Reynolds  number  3.93  x  10^ 

The  unsteady  computations  were  initialized  with  a  steady-state  solution  computed  at  the 
maximum  flapping  amplitude  where  the  flapping  velocity  is  zero,  Tuncer  and  Platzer  [Ref 
3].  See  a  copy  of  the  input  file  in  the  Appendix  A.  The  unsteady  computations  were 
carried  out  for  up  to  three  periods  of  harmonic  flapping  motion.  The  flowfields  presented 
were  taken  firom  the  last  period  of  the  computation  cycle.  All  unsteady  runs  were 
completed  at  the  Naval  Postgraduate  School  Cray  Y-MP  J94.  The  average  computer 
processing  time  (CPU)  for  15,000  time  step  of  unsteady  computations  was  9.21  hours. 

1.  Case  I 

The  flowfield  for  this  case  was  initially  computed  with  a  flapping  motion  of  A= 

.  10c  and  k=  0.5.  The  computed  streamlmes  for  both  the  steady  and  unsteady  case  are 
shown  in  Figure  12.  It  can  be  noted  in  this  figure  that  the  fluid  particles  no  longer  travel 
fi-eely  through  the  gap  to  the  upper  flap  section.  Instead,  the  suction  created  by  the 
flapping  motion  keeps  the  dividing  streamline  closer  to  the  flap  leading  edge.  The  flow 
fields  for  both  the  steady-state  and  unsteady  solution  are  shown  in  Figure  13. 
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Figure  12  Mass-Flux  Steady-State  (top)  and  Unsteady  (hot) 


Figure  13  clearly  shows  that  the  flow  is  sucked  toward  the  flap  if  the  flap  is  oscillating. 


Figure  13  Velocity  Magnitude  Computed  for  Steady-State  (top)  and  Unsteady  (bottom), 

A^O.lOc,  k=^.5 

Due  to  the  gap  size  the  effect  on  pressure  distribution  over  the  main  airfoil  was  very  small 
as  shown  in  Figure  14.  On  the  other  hand,  the  boundary  layer  profile  changes,  shown  in 
Figure  15,  at  the  aft  upper  surface  of  the  main  airfoil  are  significant.  It  is  noted  that  the 
flow  reversal  is  reduced ,  the  magnitude  of  the  velocity  vectors  is  increased  and  the 
boundary  layer  reattaches  to  the  surface  close  to  the  trailing  edge.  The  changes  are  really 
noticeable  in  the  treuling  edge  region  which  indicates  that  the  suction  produced  by  the 
flapping  airfoil  is  not  sufficient  for  a  reattachment  of  the  boundary  layer  forward  of  that 
region. 
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Figure  14  Main  Airfoil  Cp  Distribution  A  =0. 10c,  k=.  5 


x/c  =  0.80  x/c  «  0.85  x/c  »  0.90  x/c  «  0.95 


Figure  15  Boundary  Layer  Profile,  Steady  State  (top).  Unsteady  (bottom), 

A  =  0.10c,  k-.5 
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The  case  I  configuration  was  tested  at  A  =  0.05  and  k  =  .5,  with  very  similar 
results,  as  shown  in  Figure  16.  This  figure  suggests  that  due  to  the  reduced  oscillating 
amplitude  the  suction  effect  was  significantly  reduced. 


Figure  16  Velocity  Magnitude  Computed for  Steady-State  (top)  and  Unsteady  (bottom), 

A=0.05c,  k=.5 


The  time  history  of  the  unsteady  aerodynamic  loads  of  the  flapping  motion  is 
shown  in  Figure  17.  Following  the  initial  transition ,  all  aerodynamic  loads  attain  a 
periodic  behavior  with  respect  to  the  flapping  motion.  For  the  lift  coefficient  the 
oscillation  created  peak  variations  of  less  than  one  percent  for  both  cases.  Similarly,  drag 
coefficient  variations  were  less  than  2.5  percent  for  A=.  Ic  and  one  percent  for  A=.05c. 
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Figure  17  Time  History  of  the  Unsteady  Lift  and  Drag  Coefficients 
Figures  18  and  19  shows  the  particle  traces  for  four  different  flap  positions. 


2. 


Case  n 


For  this  case  the  gap  between  airfoils  was  changed  to  .25c.  Solutions  were 
obtained  for  A  =  .05  and  k  =  .5  and  1.0.  Figure  20  shows  how  changes  in  oscillating 
frequency  influence  the  thrust  generation.  As  frequency  increases,  the  suction  increases 
pulling  the  wake  of  the  main  airfoil  closer  to  the  flap  boundary. 


Figure  20  Velocity  Magnitude  Computed for  Steady-State  (top),  k  =  .5  (center)  and 

k=1.0  (bottom) 
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rv.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

Subsonic  flow  over  stationary  airfoU/oscillating  flap  combinations  were  computed 
using  a  Navier-Stokes  solver.  The  computed  flowfields  confirm  the  observations  made  in 
recent  experiments  performed  by  Dohring  [Ref.  4].  However  additional  computations  for 
other  configurations  at  different  fi-equencies,  amplitudes,  angle-of-attack,  airfoil  sections, 
etc.,  are  required  before  more  quantitative  conclusion  can  be  drawn.  Preliminary  review 
of  the  computed  boundary  layer  characteristics  indicates  that  with  a  proper  configuration 
and  applying  a  proper  combination  of  the  flapping  parameters  mentioned  above  a 
reattachment  of  the  boundary  layer  may  result. 

B.  RECOMMENDATIONS 

To  folly  satisfy  the  goals  of  this  research  a  wider  parametric  study  is  required.  In 
order  to  obtain  a  clear  picture  of  the  flapping  effects  and  to  be  able  to  quantify  the  thrust 
eflBciency  the  effect  of  the  following  parameters  hsted  below  needs  to  be  identified. 

•  Frequency  of  oscillation 

•  Amplitude 

•  Reynolds  number 

•  Gap  size  and  flap  location 

•  Angle-of-attack 

•  Airfoil  section  (i.e.  symmetric,  non-symmetiic,  flat  plate,  etc.  ) 

In  addition; 

1 .  The  complexity  of  the  flow,  especially  in  the  wake  area  of  both  airfoils, 
requires  special  attention.  For  complex  turbulent  flows,  the  choice  of  turbulence  model 
will  have  a  significant  impact  on  accuracy  and  computational  time.  Therefore,  other 
turbulence  model,  such  as  k-s  or  Baldwin-Barth,  needs  to  be  investigated  for  the  multi¬ 
element  configuration.  Nelson  [Ref  9]. 
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2.  Grid  generation  is  critical  toward  obtaining  accurate  solutions.  Tools  currently 
available  are  complex.  Therefore,  more  work  is  required  to  create  a  grid  that  solves  the 
overhang  and  overlap  problem  discussed  in  the  grid  generation  section.  A  multi-block 
grid  may  solve  this  problem. 
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APPENDIX  A 


A.  DESCRIPTION  OF  NS  INPUT  FILE  PARAMETERS 


IREAD 


ITER- 

NPRINT 

NLOAD 

ODVAR 

ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 

PLUNGE 

PLMX 

PLMY 


0  No  initial  solution,  free  stream  conditions  initialize  the  flowfield  in 
the  startup  steady-state  computations 

1  Initial  solution  is  read  from  a  binary  file  saved  from  the  previous  run 
(default,  at  the  end  of  each  run,  solution  file  is  saved  as  binary) 

2  Initial  solution  is  read  as  formatted  (plotSd  form) 

-1  Initial  solution  is  read  as  binary,  unsteady  motion  starts 

-2  Initial  solution  is  read  as  formatted  (plotSd  form),  unsteady  motion 

starts 

#  of  timesteps 

Residuals  are  printed  out  at  every  nprint  timesteps 
Aerodynamic  loads  are  printed  out  at  every  nload  timesteps 

Solution  variables,  q  array,  are  written  into  "qp.d"  file  at  every  delta  odvar 
change  in  unsteady  motions  (degrees  in  oscillatory  motion,  amplitude 
change  in  plunge) 

Steady  state  AO  A  (do  not  set  it  to  zero,  instead  set  it  to  0.0001) 
false  N/A 

true  sinusoidal  oscillations  in  pitch 
false  N/A 

true  straight  ramp  motion  in  pitch 

Reduced  frequency  of  the  unsteady  pitching  motion  based  on  the  half 
chord,  chord  length  is  assumed  to  be  1. 

Min  AO  A  of  the  pitching  motion 

Max  AOA  of  the  pitching  motion 

false  N/A 

Plunge  amplitude  in  x 
Plunge  amplitude  in  y 
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PLPHSXD 

PLFREQ 

REYNOLD 

MACH 

Vise 

TURBL 

TRANS 


TIMEACC 


COUR 


NEWTIT 

POTEN 


The  Phase  angle  between  x  and  y  amplitudes,  in  degrees 

The  reduced  frequency  of  the  plunging  motion 
based  on  the  half  chord,  chord  length  is  assumed  to  be  1 . 

Reynolds  number  of  the  freestream  flowfield 

Mach  number  of  the  freestream  flowfield 

false  Euler  solution 

true  Viscous  Navier-Stokes  solution 

false  Laminar  flow  is  assumed. 

true  Baldwin-Lomax  turbulence  model  is  applied. 

false  Fully  turbulent  flow  is  assumed. 

true  BL  transition  with  Michel  criterion  is  modeled.  A  steady-state 
solution  has  to  be  obtained  first.  TIMEACC  has  to  be  set  to  true. 

false  Variable  local  time  stepping  in  the  computational  grid,  used  in  the 
computations  of  steady  state  and/or  attached  flowfields. 

true  Constant  time  stepping  everywhere  in  the  computational  grid,  used 
in  the  computations  of  unsteady  and  separated  flowfields 

Courant  number  of  the  timestepping  (50-1500),  determines  the  timestep  of 
the  computations  based  on  the  minimum  grid  size  and  the  freestream 
conditions.  Its  value  depends  on  the  computational  grid.  For  diverging 
computations,  when  the  residuals  in  the  output  file  increases  in  time,  it  is 
the  sign  that  its  value  is  to  be  reduced. 

Number  of  Newton  subiterations  in  each  timestep,  applied  in  unsteady 
flows  (2-3),  for  steady  flowfields  it  is  set  to  1 . 

This  variable  is  reserved  for  NS-Potential  flow  interactive  solution 
procedure. 


The  NS  program  needs,  in  general,  the  following  input  files: 


ns.in:  Input  file  which  defines  certain  parameters  as  given  above.  It  is  set  up  for  starting 

a  generic  steady-state  solution.  For  continuing  computations  or  for  any  changes 
in  the  input  variables  you  need  to  edit  this  file. 

grid.in  :  A  formatted  file  in  which  the  computational  grid  coordinates  are  stored. 

strs.d  ;  A  file  in  which  the  starting  solution  is  stored.  Could  be  binary  or  formatted  in 
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accordance  with  the  IREAD  variable. 

Once  you  have  the  input  files  ready,  you  can  run  the  NS  solver  on  SGI's  by  simply 
submitting  it  with  the  following  command  at  the  prompt; 

run  ns 


B.  STEADY-STATE  CONDITIONS 

C..  IREAD,  NITER  NPRINT,  NLOAD  ODVAR 
0  3000  10  10  0.01 

C..  ALPHA  OSCIL  RAMP  REDFRE  ALFAMND  ALFAMXD 
10.0  false  false  0.099  0.001  10.0 

C..  PLUNGE  PLMX  PLMY  PLPHSXD  PLFREQ 
false  0.  0.10  0.  0.5 

C..MACH  REYNOLD  VISC  ITURBL  TRANS 

0.300  3930000.  true  1  false 

C..  TIMEACC  COUR  NEWTIT 
false  100.  1 

C.  UNSTEADY  CONDITIONS 

C..  IREAD,  NITER  NPRINT,  NLOAD  ODVAR 
-2  13000  20  20  0.01 

C..  ALPHA  OSCIL  RAMP  REDFRE  ALFAMND  ALFAMXD 
10.0  false  false  0.099  0.001  10.0 

C..  PLUNGE  PLMX  PLMY  PLPHSXD  PLFREQ 
true  0.0  0.10  0.0  0.5 

C..MACH  REYNOLD  VISC  ITURBL  TRANS 

0.300  3930000.  true  1  false 

C..  TIMEACC  COUR  NEWTIT 
true  500.  1 


D.  GRMOVE  (MOTION  SUBROUTINE) 


subroutine  grmove(dalfa,dx,dz) 
include  'corns,  f 

dimension  xold(nia,nka),  zold(nia,nka) 

do  i=l,imx(l) 

do  k=l,kmx(l) 

xold(i,k)  =  x(i,k) 

zold(i,k)  =  z(i,k) 

enddo 
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enddo 


if(  dalfa  .ne.  0.)  then 
ca  =  cos(  dalfa ) 
sa  =-sin(  dalfa  ) 
do  i=l,inix(l) 
do  k=l,kmx(l) 

x(i,k)  =  xold(i,k)  *  ca  -  zold(i,k)  *  sa 
z(i,k)  =  zold(i,k)  *  ca  +  xold(i,k)  *  sa 
enddo 
enddo 
endif 

IF(  dx  .ne.  0.  or.  dz  ne.  0. )  then 
if(.not.  distort)  then 
do  k=l,knix(l) 
do  i=l,imx(l) 
x(i,k)  =  x(i,k)  +  dx 
z(i,k)  =  z(i,k)  +  dz 
enddo 
enddo 
else 

kiimer  =  45 
kouter  =  80 
iinnerl  =89 
iinner2  =283 
iouterl  =108 
iouter2  =264 
xrangeinner  =  3.5 
xrangeouter  =  2. 

yrangeinner  =  2. 
yrangeouter  =  2. 
do  i=l,imx(l) 

if  (i  .le.  iinnerl  .and.  i  .ge.  1)  then 
xfactor  =  1 

elseif  (i  .ge.  iinner2  and.  i  le.  inix(l))  then 
xfactor  =  1 

elseif  (i  .It.  iouterl  .and.  i  .gt.  iinnerl)  then 
slope  =  -(yrangeinner+yrangeouter)/ 

>  float((iouterl-l)-(iinnerl+l)) 

xcept  =  yrangeinner-slope*float(iinnerl+l) 
yinput  =  slope*float(i)+xcept 
xfactor  =  0.5*(tanh(yinput)+l .) 
elseif  (i  gt.  iouter2  and.  i  .It.  iinner2)  then 
slope  =  (yrangeinner+yrangeouter)/ 

>  float((iinner2-l)-(iouter2+l)) 
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xcept  =  yrangeinner-slope*£loat(iinner2+l) 
yinput  =  slope*float(i+2)+xcept 
xfactor  =  0.5*(tanh(yinput)+l .) 
else 

xfactor  =  0. 
endif 

do  k=l,lanx(l) 
k  .le.  kinner )  then 
factor  =  1. 

elseif(  k  .gt.  kinner  .and.  k  .It.  kouter )  then 
slope  =  -(xrangeinner+xrangeouter)/ 

>  float((kouter-l)-(kinner+l)) 

xcept  =  xrangeinner-slope*float(kinner+l) 
xinput  ==  slope*float(k)+xcept 
factor  =  0.5*(tanh(xinput)+l.) 
else 

factor  =  0. 
endif 

xold(i,k)=  x(i,k) 
zold(i,k)=  z(i,k) 

x(i,k)  =  x(i,k)  +  dx 
z(i,k)  =  z(i,k)  +  dz*factor  *  xfactor 
enddo 
enddo 
endif 
ENDBF 

ifl^  unstdy )  then 

do  k=l,knix(l) 
do  i=l,inix(l) 

xtau(i,k)  =  (x(i,k)  -  xold(i,k))/dtau 
ztau(i,k)  =  (z(i,k)  -  zold(i,k))/dtau 
enddo 
enddo 
else 

do  k=l,knix(l) 
do  i=l,imx(l) 
xtau(i,k)  =0. 
ztau(i,k)  =0. 
enddo 
enddo 
endif 

call  metric 
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if(  not.  unstdy  or.  itr  .eq.  niter )  then 
open(unit=9 1  ,file='grid .  out',form- formatted') 
write  (91,*)  imx(l),  kmx(l),  iwks(l),  iwall 
write  (91,*)  ((x(i,k),  i  =  1,  imx(l)),  k  =  l,kinx(l) ), 
>  ((z(i,k),  i  =  1,  imx(l)),  k  =  l,kmx(l) ) 

endif 

return 

end 
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APPENDIX  B 


A.  GRIDGEN  INPUT  FILE  (NACA  4412) 


100 

1.833333 

0 

0 

1.832997 

-1.4E-05 

0 

1.831988 

-5.7E-05 

0 

1.830313 

-0.00013 

0 

1.827978 

-0.00023 

0 

1.824994 

-0.00036 

0 

1.821375 

-0.00051 

0 

1.817136 

-0.0007 

0 

1.812297 

-0.00091 

0 

1.806878 

-0.00116 

0 

1.800905 

-0.00143 

0 

1.794402 

-0.00173 

0 

1.787397 

-0.00206 

0 

1.779921 

-0.00243 

0 

1.772005 

-0.00282 

0 

1.763682 

-0.00325 

0 

1.754987 

-0.00371 

0 

1.745956 

-0.00421 

0 

1.736626 

-0.00473 

0 

1.727034 

-0.00528 

0 

1.717219 

-0.00586 

0 

1.707221 

-0.00646 

0 

1.697078 

-0.00708 

0 

1.68683 

-0.00771 

0 

1.676518 

-0.00834 

0 

1.66618 

-0.00896 

0 

1.655855 

-0.00956 

0 

1.645584 

-0.01013 

0 

1.635403 

-0.01066 

0 

1.625516 

-0.01116 

0 

1.615838 

-0.01167 

0 

1.60636 

-0.01216 

0 

1.597109 

-0.01264 

0 

1.588117 

-0.01307 

0 

1.579412 

-0.01343 

0 

1.571021 

-0.01371 

0 

1.562971 

-0.0139 

0 

1.55529 

-0.01397 

0 

1.548004 

-0.01391 

0 

1.541138 

-0.01371 

0 

1.534718 

-0.01336 

0 

1.528769 

-0.01284 

0 

1.523314 

-0.01216 

0 

1.518378 

-0.pl  13 

0 

1.513983 

-0.01025 

0 
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1.510148 

-0.00902 

0 

1.506894 

-0.0076 

0 

1.504238 

-0.00599 

0 

1.502195 

-0.00419 

0 

1.500779 

-0.00219 

0 

1.5 

0 

0 

1.500433 

0.00471 

0 

1.501666 

0.007159 

0 

1.503578 

0.009655 

0 

1.506166 

0.012185 

0 

1.509425 

0.014728 

0 

1.513346 

0.017264 

0 

1.517917 

0.019767 

0 

1.523122 

0.022211 

0 

1.528943 

0.02457 

0 

1.535358 

0.026815 

0 

1.54234 

0.028917 

0 

1.549861 

0.030851 

0 

1.557887 

0.032591 

0 

1.566384 

0.034114 

0 

1.575313 

0.035397 

0 

1.584631 

0.036424 

0 

1.594297 

0.03718 

0 

1.604266 

0.037653 

0 

1.614489 

0.037837 

0 

1.624921 

0.037726 

0 

1.63547 

0.037324 

0 

1.645972 

0.036695 

0 

1.656548 

0.035879 

0 

1.667154 

0.034885 

0 

1.677746 

0.033726 

0 

1.688281 

0.032416 

0 

1.698716 

0.03097 

0 

1.709009 

0.029403 

0 

1.71912 

0.027731 

0 

1.729007 

0.025973 

0 

1 .738634 

0.024145 

0 

1.747962 

0.022265 

0 

1.756955 

0.020352 

0 

1.76558 

0.018425 

0 

1.773803 

0.016503 

0 

1.781595 

0.014604 

0 

1.788926 

0.012748 

0 

1.795769 

0.010954 

0 

1.802101 

0.009241 

0 

1.807898 

0.007626 

0 

1.813139 

0.006128 

0 

1.817806 

0.004762 

0 

1.821884 

0.003544 

0 

1.825358 

0.002489 

0 

1.828216 

0.001607 

0 
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1.830449  0.00091  0 

1.83205  0.000407  0 

1.833012  0.000102  0 

1.833333  0  0 

B.  ELLIPTICAL  SOLVER  PORCESS 

•  Select  Domain  commands  from  GRIDGEN’s  main  menu. 

♦  Select  the  domain  to  be  solved. 

•  Set  solver  attributes  (Boundary  conditions) 

♦  Foreground  control  function 

=>  Enable  surface  edges  (ON) 

=>  Edge  influence  (ON) 

Enter  AS  delay  factor  =10.  (Number  of  gridlines 

from  the  surface  for  which  these  conditions  apply.) 
=?>  Enter  angle  decay  factor  =10.  (Angle  gridlines 
makes  with  the  boundary  itself) 

♦  Background  control  function 

Select  Laplace  function 

♦  Done  with  settings 

♦  Run  solver. 
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INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center . 2 

8725  John  J.  Kingman  Rd.  STE  0944 

Ft.  Belvoir,  VA  22060-6218 

2.  Dudley  Knox  Library  . . 2 

Naval  Postgraduate  School 

4111  Dyer  Rd. 

Monterey,  CA  93943-5101 

3.  Mr.  MaxF.  Platzer 

Department  of  Aeronautics  and  Astronautics,  Code  AA/Pl . 10 

Naval  Postgraduate  School 
Monterey,  CA  93943-5101 

4.  Mr.  Ismail  H.  Tuncer 

Department  of  Aeronautics  and  Astronautics,  Code  AA/Tu  . 1 

Naval  Postgraduate  School 
Monterey,  CA  93943-5101 

5.  LCDR  Miguel  A.  Ortiz . 2 

105  ShubrickRoad 

Monterey,  CA  93940 
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